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SECTION  I 


INTEGRATED  AEROSPACE  ENGINE  MANAGEMENT 

1.  DESCRIPTION  OF  THIS  REPORT 

Effective  management  of  high  cost  equipment  requires 
accurate  information  about  the  operating  lifetimes  of  the 
equipment.  There  are  few  items  except  aerospace  engines  for 
which  it  is  worthwhile  to  collect  data  on  operating  or  flying 
time.  Since  the  trouble  is  taken,  the  best  possible  infor- 
mation should  be  derived  from  that  data,  and  that  information 
should  be  used  to  the  fullest  extent  possible  within  cost  and 
procedural  limitations. 

The  potential  for  an  integrated  system  for  aerospace 
engine  planning,  managing,  and  decision  making  may  be  suf- 
ficient to  deserve  consideration  for  some  changes  in  procedures. 
The  analytical  and  statistical  aspects  of  the  current  pro- 
cedure, the  actuarial  technique,  and  the  potential  for  a 
future  engine  management  system  are  described  in  this  section. 
The  remainder  of  the  report  is  to  provide  some  foundations 
I for  the  future  system,  foundations  in  estimation  and  prediction 

of  engine  removals. 

The  work  statement  leading  to  this  report  essentially 
requested  more  information  from  less  data.  That  is  a chal- 
lenge to  any  statistician,  and  often  the  way  to  improve  in- 
formation output  is  to  incorporate  newer  techniques. 
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Unfortunately,  these  advancements  may  be  accompanied  by  in- 
creased data  and  computation  requirements.  Not  so  in  the 
circumstances  of  this  report.  In  Section  II,  a maximum  like- 
lihood estimator  of  the  distribution  function  of  flying  hours 
till  engine  removal  is  derived.  This  estimator  may  yield 
more  information  from  less  data  than  is  now  needed  for  the 
actuarial  method,  and  for  less  computational  effort.  In 
Section  III  is  a plan  which  may  yield  still  more  accuracy. 

Also  fundamental  to  an  integrated  engine  management 
system  is  a prediction  of  engine  replacement  requirements. 

The  simulation  procedure  for  prediction  of  replacements  may 
be  used  with  the  actuarial  estimator  or  with  the  estimator 
obtained  in  Section  II  or  III. 

Section  V is  an  annotated  bibliography  chosen  for 
relevance  to  engine  removal  time  estimation  and  to  integrated 
engine  management  systems  in  the  future. 
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2.  THE  ACTUARIAL  TECHNIQUE 

Thtt  statistical  procedure  now  used  in  aerospace  engine 
management  is  based  on  the  actuarial  technique  originally 
developed  for  use  by  life  insurance  companies.  At  the  time 
of  its  adoption  by  the  Air  Force,  there  was  little  choice  of 
method  for  estimation  of  engine  lives,  and  an  admirable  job 
was  done  of  adapting  the  technique  to  Air  Force  needs.  The 
actuarial  technique  is  the  standard  for  comparison  of  esti- 
mation techniques,  and  any  proposed  technique  must  provide  the 
same  form  of  information  from  the  seune  data  so  there  is  as  little 
disruption  as  possible  if  there  is  eventual  change.  Therefore, 
it  is  appropriate  to  describe  the  actuarial  technique  and  its 
inherent  limitations. 

The  actuarial  technique  as  applied  to  aerospace  engine  sys- 
tems is  used  to  obtain  failure  rate  information  from  data  on 
flying  hours  till  engine  replacement.  This  is  in  addition  to 
accounting  and  inventory  data  on  engines.  The  statistical 
part  of  the  information  is  svjnmarized  in  a quarterly  report 
to  the  Aerospace  Engine  Life  Committee.  This  Committee 
monitors  the  performance  of  each  engine  type  and  its  main- 
tenance system  of  base  and  depot  repair  shops  and  inventories. 
Extensive  reports  are  provided  quarterly  to  the  Air  Force 
Logistics  Command,  and  the  information  may  also  be  distributed 
and  used  in  ways  unknown  to  the  author. 


Data  on  each  engine  is  in  the  form  of  flying  hours  till 
engine  replacement  and  the  replacement  cause.  There  is  also 
a quarterly  report  on  accumulated  flying  hours  of  engines  in 
service  and  other  information  of  an  accounting  nature.  Let  T 
denote  the  random  variable,  flying  hours  till  replacement,  for 
a particular  type  of  engine  and  a specified  set  of  removal  codes. 

The  so-called  "removal  rates"  estimated  by  the  actuarial  method 
are  the  probabilities  of  removal  in  a flying  hour  interval 
given  the  engine  has  survived  till  that  interval , 

P{Te [kAt, (k+1) At] lT>kAt} , k=0 , 1 , 2 , . . . ,m  (1) 

where  At  is  a conveniently  chosen  interval  of  from  1/30  to  1/100 
of  the  maximum  allowable  flying  time  of  an  engine.  The  choice  of 
At  determines  m.  A simple  estimator  is  the  ratio  of  the  number 

of  engines  that  are  replaced  in  an  interval  to  the  number  that 
survived  replacement  to  the  beginning  of  that  interval.  The 
actuarial  technique  also  incorporates  information  about  engines 
that  may  have  been  flown  partway  through  an  interval  at  the 
time  of  the  quarterly  report. 

The  actuarial  technique  was  originally  developed  for  esti- 
mation of  human  lifetimes,  so  rather  large  fluctuations  in  the 
estimator  of  (1)  were  unsatisfying  in  appearance.  There  is 
no  inherent  reason  why  the  failure  rate  of  human  beings  should 
change  drastically  from  one  interval  to  the  next,  therefore  the 
ratio  estimators  of  (1)  are  smoothed  by  the  moving  average 

method.  Where  dat^  is  scarce,  the  smoothed  estimator  is  extrapolated 
to  the  maximum  flying  time  by  linear  regression.  The  resulting 


T. 


estimator  of  the  sequence  of  Interval  replacement  rates  (1)  Is 


graphed  by  hand  and  presented  to  the  Aerospace  Engine  Life 


Committee  along  with  several  other  measures  of  performance  of 


the  engines  and  the  repair  system. 


The  actuarial  technique  as  applied  in  the  Air  Force  is 


remarkably  well  conceived,  planned,  and  executed  considering 


the  era  in  which  the  procedure  was  adopted,  the  late  1950 's. 


Since  then,  some  of  the  calculation  has  been  progreunmed  for 


computers,  and  automated  data  collection  may  be  implemented 


for  some  engine  types.  But  it  appears  to  the  author  that 


the  procedure  and  the  use  of  the  engine  replacement  informa- 


tion has  not  varied  significantly  since  the  procedure  was 


first  used. 


There  is  only  one  fundamental  criticism  of  the  actuarial 


technique  as  applied  to  aerospace  engines.  The  graphs  of  the 


replacement  rates  are  far  from  smooth.  Sharp  peaks  appear  in 


the  neighborhood  of  periodic  inspection  times  despite  the 


moving  average  smoothing.  There  is  information  in  the  fact 


that  many  engine  replacements  occur  at  the  times  of  periodic 


inspection,  and  this  information  is  exploited  in  the  estimator 


derived  in  Section  II. 


Other  events  since  the  implementation  of  the  actuarial 


technique  increase  the  pressure  for  development  of  new 


methods  to  get  more  information  from  less  data.  Advances  in 


statistics  and  reliability  theory  have  occurred.  Prolifera- 


tion of  engine  types  has  resulted  in  less  data  for  each  type. 


f- 


hampering  the  use  of  the  actuarial  technique  which  depends 
on  having  a large  amount  of  data.  The  future  engine  manage - 
ment  decisions  about  logistics  given  a required  flying  hour 
program  may  require  more  accurate  and  different  information 
about  engine  performance.  The  data  exists  now;  it  must 
be  exploited  more  fully#  but  without  disruption  of  existing 
procedures . 


i 


i 

\ 
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3.  PROPOSAL  FOR  AN  INTEGRATED  ENGINE  MANAGEMENT  SYSTEM 
That  part  of  the  Air  Force  associated  with  aerospace 
engines  may  be  considered  as  a multi-product,  multi-echelon 
production  system.  The  products  are  flying  hours  for  each 
engine  type.  This  analogy  allows  the  application  of  analyti- 
cal decision  making  techniques  developed  for  other  organiza- 
tions. The  methods  of  estimation  and  prediction  developed 
in  Sections  II,  III  and  IV  provide  bette’'  information  funda- 
mental to  decision  making.  That  information  can  penetrate 
the  engine  management  system  along  existing  paths  to  the 
point  where  it  is  needed  for  analytical  decision  making. 
Parallel  in  time  to  the  penetration  of  better  information 
and  the  analysis  for  decision  making,  a major  effort  should 
be  made  to  integrate  the  decision  making  policies  throughout 
the  Air  Force,  vertically  in  the  engine  management  system 
and  horizontally  across  activities  such  as  civilian  and 
military  personnel  management,  missile  system  management, 
etc.  What  follows  is  a description  of  the  engine  management 
system  considered  as  an  integrated  production  system. 


] 

I 


1 


1 

1 

il 


Aerospace  engine  management  may  be  analyzed  as  management 
of  a production  system.  The  product  is  flying  hours.  The 
sales  manager,  HQ  USAF,  predicts  the  required  number  of  fJvina 
hours  for  each  command  and  this  becomes  a production  require- 
ment which  is  further  broken  down  into  engine  type 
flight  requirements.  Headquarters  USAF  and  the  entire 
logistics  system  necessary  to  support  the  engines  constitute 
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the  engine  management  system  referred  to  throughout  this 
report. 

H.  M.  Wagner  [1]  describes  an  integrated  approach  to  the 
management  of  multi-product,  multi-echelon  production  systems. 

He  focuses  on  nine  aspects  of  the  production  system  that  must 
be  considered  if  change  in  any  one  is  contemplated.  Translated 
into  the  context  of  the  aerospace  engine  management  system,  these 
foci  are: 

1)  forecasting  demand  for  and  flight  capabilities  of 
the  entire  system  throughout  the  planning  horizon, 

2)  allocation  of  the  flying  hour  program  among  commands, 
bases , and  engine  types , 

3)  purchasing  and  overhaul  of  engines  and  spare  ptrts, 

4)  inventory  levels  at  every  point  within  the  system, 

5)  future  target  inventory  levels  to  cover  contingencies 
and  anticipated  flight  requirements, 

6)  backlog  policies  to  handle  shortfalls  in  flying 
programs  and  inventories, 

7)  the  management  information  system, 

8)  the  contingencies  in  the  flying  programs,  and 

9)  the  management  system,  its  structure,  policies, 
and  measures  of  effectiveness. 

At  each  of  these  foci  there  are  decisions  that  may  re- 
quire information  about  engine  removal  times  and  the  number 
of  replacements  required  to  achieve  a specified  flight  plan. 
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This  information  should  penetrate  the  logistics  system  to 
the  point  needed  and  in  the  form  needed.  Analysis  for 
decision  making  will  help  determine  the  point  of  need  and 
the  form  of  the  information  required  for  decisions.  Even- 
tual integration  will  take  place  as  relevant  information 
infiltrates  the  engine  management  system  and  as  analysis 
shows  the  relationship  of  decisions  made  anywhere  in  the 
engine  management  system. 

There  are  many  ways  that  information  about  engine 
removal  times  and  prediction  of  replacements  enters  into 
analysis  for  decision  making  at  all  nine  foci. 

1)  Clearly,  a prediction  of  the  number  of  replacements 
necessary  to  achieve  a flying  program  is  relevant  to 
forecasting  the  capability  to  meet  planned  flight 
requirements.  Also,  estimates  of  the  engine  re- 
moval times  are  relevant  because,  through  inter- 
actions between  the  engine  age  in  its  product  life 
cycle  and  maintenance,  inventory,  and  purchasing 
policies , it  is  possible  to  be  left  without  enough 
engines  to  meet  contingencies.  Risk  analysis  should 
determine  the  probability  of  meeting  contingencies 
throughout  the  product  life  cycle  of  an  engine  type. 

2)  Allocation  of  the  flying  program  follows  from  the  fore- 
cast of  requirements  and  capabilities.  Mathematical 
programming  will  assist  allocation,  but  many  other 
factors  must  be  recognized  before  final  allocation 

is  made. 
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3)  Purchase  of  engines  and  spare  parts  and  overhaul 
decisions  depend  on  the  future  flight  requirements 
and  the  performance  of  the  engines  as  Indlcatrd  by 
the  prediction  of  the  number  of  replacements  to  meet 
the  requirements. 

4)  Optimal  engine  inventory  policies  may  be  determined 
analytically  to  control  stockouts  and  costs . Parts 
inventory  policies  follow  from  engine  policies.  The 
demand  for  engines  as  predicted  by  the  simulation 
model  in  Section  IV  or  other  means  is  a necessary 
input  for  inventory  policy  analysis. 

5)  Future  target  inventory  levels  are  determined  indirectly 
from  anticipated  flying  programs,  contingencies,  and  the 
capabilities  of  the  engines  themselves.  Some  redistri- 
bution may  be  necessary  to  meet  the  targets,  and  again, 
mathematical  programming  techniques  can  help  analysis 
for  redistribution. 

6)  Alternate  policies  for  making  up  backorders  in  inven- 
tories and  shortfalls  in  meeting  the  flying  hour  pro- 
gram to  determine  their  effect  on  other  parts  of  the 
engine  management  system. 

7)  The  management  information  system  should  provide  in- 
formation in  the  form  needed  at  the  point  necessary 
for  decision  making  in  foci  1-6.  An  information  sys- 
tem for  engine  replacements  already  exists;  additional 
or  more  accurate  information  may  be  necessary  at  some 
points  but  not  at  other  points.  The  information 
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gathering  system  should  also  be  considered  to  deter- 
mine whether  the  effort  expended  is  worthwhile. 

And  work  such  as  that  represented  by  this  report  should 
be  done  to  extract  more  information  from  less  data. 

8)  Contingencies  in  the  flying  progrcun  may  result  from 
external  or  internal  causes.  Internal  contingencies 
may  be  predicted  by  risk  analysis  of  policies  in  foci 
1-6. 

9)  The  military  command  system  provides  a management 
structure  unmatched  by  any  other  production  system. 

This  structure  has  inherent  inertia  that  calls  for  more 
justification  for  change  than  in  a business  where  the 
consequences  of  risk  are  not  so  large.  Any  changes 
should  be  accomplished  in  parallel  with  existing  pro- 
cedures. The  estimation  procedure  in  Section  II  is 
designed  for  use  along  with  the  actuarial  method 
until  parts  of  the  actuarial  method  may  be  replaced. 
Analysis  for  decision  making  may  proceed  at  all  levels 
of  the  engine  management  system  and,  if  results  are 
worthwhile,  information  necessary  for  the  analytical 
decision  making  should  be  made  available.  Meanwhile, 
integration  of  decision  making  throughout  the  engine 
management  system  is  facilitated  by  the  command 

and  communications  systems. 


This  completes  a description  of  the  context  and  anti- 
cipated development  of  integrated  engine  management.  The 
contributions  in  the  remaining  sections  are  foundations  for 
obtaining  better  information  for  use  in  analytical  decision 


making . 


SECTION  II 

MAXIMUM  LIKELIHOOD  ESTIMATION  OP  THE 
REMOVAL  TIME  CUMULATIVE  DISTRIAUTIOM  FUNCTION 


1 . INTRODUCTION 

The  problem  is  to  estimate  the  distribution  function  of 
engine  removal  times  which  is  assumed  to  have  the  form 

F(t)  - 1 - r,(t)  • n (1-p.)  (2) 

^ i«l  ^ 

for  t ^ 0.  An  estimate  is  necessary  to  predict  demands  for  re- 
placement engines.  The  form  of  the  distribution  function  is 
suggested  because  of  the  high  proportion  of  removals  that  occur 
at  inspections.  Each  p^  is  the  probability  of  remo.al  at  the 
time  of  the  i^^  inspecfion,  i»l,2,...,  k£».  The  limit  i^^j 
is  the  index  of  the  last  inspection  prior  to  or  at  time  t. 

The  product  of  1-p^  terms  is  assumed  equal  to  1 if  there 
are  no  inspections  at  or  before  t.  The  tail  distribution 
Fj^(t)  * l-Fj^(t)  is  the  probability  distribution  of  usage 
removal  time.  That  distribution  may  be  truncated  at  a maxi- 
mum operating  time,  when  ail  engines  are  removed  for 

overhaul.  The  model  (2)  is  also  appropriate  to  the  life  test- 
ing situation  with  all  items  put  on  test  at  one  time  and  a 
random  number  of  survivors  removed  from  test  at  known  times 
prior  to  termination  of  the  test  (progressive  censoring) . 

Sample  data  tell  the  engine  removal  time  and  whether 
removal  occurs  during  usage,  at  an  inspection,  or  at  the 
maximum  time.  (Actually  inspection  and  maximum  time  removals 


vary  slightly  around  the  scheduled  times,  but  fixed  inspection 
tiroes  and  t will  be  assuroed.  Estimation  of  inspection 
removal  tiroes  constitutes  another  problem.)  The  sample  data 
may  be  used  to  construct  estimators  of  F(t)  in  several  ways. 

The  empirical  distribution  function  with  jumps  inversely 
proportional  to  sample  size  occurring  at  all  removal  times  may 
be  used  to  estimate  F(t) . The  sample  may  be  segregated  into 
usage,  inspection  and  maximum  time  removals;  the  usage  removal 
times  may  be  used  to  construct  an  empirical  distribution  function 
for  Fj^(t),  and  inspection  and  maximum  time  probabilities  ir.ay 
be  estimated  from  sample  proportions. 

The  method  of  maximum  likelihood  is  used  to  determine  an 

A 

estimator  of  F(t) , denoted  F(t) , which  incorporates  the  infor- 
mation that  an  engine  removed  at  an  inspection  or  at  has 

survived  previous  usage  and  vice-versa.  That  estimator  for 
F(t)  is  the  empirical  distribution  function,  but  the  estimators 
of  the  p^,  i=l,2,...,k  and  of  Fj^Ct)  are  not  the  sample  pro- 
portions and  the  empirical  distribution  function! 

The  maximum  likelihood  method  insures  that  estimators  will 
be  functions  of  the  sufficient  statistics.  Asymptotic  pro- 
perties of  unbiasedness,  efficiency,  consistency  and  normality 
can  probably  be  proved  by  verification  of  regularity  condi- 
tions for  the  likelihood  function.  The  sample  proportions  of 
survivors  which  are  derived  to  estimate  the 
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i=l,2...,k,  may  also  have  good  small  sample  properties. 

The  Glivenko-Cantelli  theorem  may  be  applied  to  verify  con- 
sistency of  the  estimator  for  Fj^(t),  and  perhaps  the 
Kolmogorov  theorem  may  be  extended  to  find  the  asymptotic 
distribution  of  the  difference  between  F^{t)  and  its  estimator. 
All  these  results  are  suggested  because  of  the  reasonable  form 
of  the  estimator  of  F(t) . They  have  not  been  proved. 

Comparison  of  F(t)  to  other  methods  of  estimation  is  in  an 

A 

incomplete  state  also.  However,  F(t)  may  prove  to  be  better 

than  the  actuarial  estimator  which  breaks  [0,  t„^„]  into 

max 

approximately  60  equal  intervals  and  estimates  F(t)  as  piece- 
wise  exponential  on  those  intervals . 
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2.  derivation  of  the  maximum  likelihood  estimator 


The  engine  removal  times  Tj^,...,Tjj  in  the  sample  are 
assumed  to  be  independent  and  identically  distributed  accord- 
ing to  F(t) , (2).  In  addition,  the  engine  removal  code  speci- 
fies whether  the  removal  was  during  usage,  at  inspection,  or 
at  maximum  time.  Define  the  following  sequences: 

{nj},  the  number  of  removals  at  the  inspection  and 
j=l»2, , . . ,k+l; 

{m^},  the  number  of  usage  removals  between  the  (j-l)st  and 

j 

the  jth  inspection;  M4=I  m . , j=l,2, . . . ,k+l; 

•^i=l  ^ 

{tj.},  the  sequence  of  usage  removal  times,  i=l,  2 , . . , ,Mj^+2.» 
{tj},  inspection  times  and 

These  sequences  turn  out  to  be  sufficient  statistics.  The 
sample  size  N can  be  expressed  as 


N 


k+1 

E 

j = l 


n 


j 


in  terms  of  usage  removals,  inspection  removals,  and  maximum 
time  removals. 
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The  likelihood  function  for  the  sample  information  is 


with  denoting  the  density  of  Fj^(t),  discrete  or  continuous, 

^max  " ^max"^  small  e>0,  and  the  awkward  notation  j is 
used  to  indicate  the  index  of  the  last  inspection  prior  to  the 
usage  removal  at  t^.  The  first  term  of  the  likelihood  function 
is  the  probability  of  the  usage  removals,  the  second  term  is  the 
probability  of  inspection  removals,  and  the  last  is  the  proba- 
bility of  all  the  survivals  to  t . C is  a combinatorial  constant. 
^ max 

In  order  to  make  the  likelihood  function  positive,  all  of 

the  fjCtp  should  be  positive,  and  to  maximize  it,  l-Fj^(tj)  and 

1-Fj^  should  be  as  small  as  possible  consistent  with  the 

constraints  on  the  distribution  function  F(t);  F(0")»0,  F(»)*l, 

and  F(t)  non-decreasing  in  t.  As  in  the  maximum  likelihood 

derivation  of  the  e-rpirical  distribution  function,  all  the  mass 

of  F,  (•)  should  be  placed  at  the  observations  {t  } and  t 

1 ^ 1 max 

Mj 

1-Fi(t  ) - 1-  I fi(tj) 


Then 


However,  not  all  of  the  jumps  in  Fj^(-)  will  be  of  the  same  size 
as  in  the  empirical  distribution  function.  The  derivative  of 
the  natural  logarithm  of  the  likelihood  function  with  respect 
to  gives 

M. 

61nL/6f(tp  « l/fl(t.)  - fj(t.)/[l- 

i=l 


for  all 


, j*l,2, . . . ,k+l.  This  indicates  that 


the  maximum  likelihood  estimators  of  the  jumps  will  be  the  same 
in  the  intervals  before  the  first  inspection,  between  inspec- 
tions, and  after  the  last  inspection.  These  jump  sizes  will  be 
denoted  by  f , f 2 , . . . , f . The  likelihood  function  now  simpli- 


fies to 
^1 


(k^ 

(x=l 


i-1 

f.  n ci-p,) 

" j=i 


j 


M .M 
i i-1 


( k+1 
' n 


p.  n ci-p.) 

3 i=i  1 


It  is  a fairly  long  algebraic  exercise  to  show  that  61nL/6f^  = 

0,  j *1 , 2 , . . . k+1 , and  61nL/6pj  = 0,  j=l,2,...,k  are  satisfied  by 


f.  . N 


-1 


and 


Pj  = "j 


j-1 

n (i-f5.) 
i=l  3 

j-J 


-1 


r 3;^  1 

N-M.-  y n. 

3 i = l 3 


-1 


Insertion  of  these  results  into  equation  (2)  shows  F(t)  to 
be  the  empirical  distribution  function'. 

Although  the  conditions  for  maxima  on  the  matrix  of  second 
derivatives  of  the  likelihood  function  have  not  yet  been 
checked,  the  reasonableness  of  the  estimators  suggests  they 
are  maximum  likelihood  rathei  than  minimum  likelihood.  These 
estimators  are  non-negative  and  satisfy  constraints  derived 
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from  the  constraints  on  P(t) , so  they  will  be  presumed  to  be 
maximum  likelihood  estimators.  The  form  of  the  estimators 

A 

has  some  intuitive  explanation.  The  estimators  p j , j=l,2,...,k, 
are  simply  the  proportions  of  inspection  removals  at  the 
inspection  out  of  all  the  engines  that  survived  to  that 

A 

inspection.  The  f ^ are  the  usual  jumps  o6  the  empirical 
distribution  function,  1/N,  but  conditioned  on  survival 
through  all  previous  inspections.  The  estimator  of  Fj^(t)  is 

A 

shown  in  Figure  1.  Notice  that  the  jumps  in  Fj^(t)  are  larger 
and  more  sparse  in  the  later  inter-inspection  intervals 
because  fewer  engines  have  survived  inspections . 


other  estimators  of  potential  interest  are  available 
or  are  easily  obtained  from  the  maximum  likelihood  estimators. 


The  probability  that  an  engine  will  be  removed 
while  in  use  is 


1 max 


k+1. 

) = 1 f.m. 

j=l  3 D 


-The  probability  of  survival  to  maximum  time  is 

1-F(t"  ) = n,  .^N 

max  K+l 

- The  failure  rate  function  r(t)  is  obtainable  from 
the  relation 

t^ 

l-F(t)  = exp(- / r(u)du) 

0 

- The  probability  of  removal  in  any  interval  of  size 


w given  survival  to  time  t is 


F (t+w) 

F(t+w)-F(t)  ^ V__  _ n 

A A 


1*-F(t) 


1=1  t 


(1-Pi) 


(These  estimates  would  replace  the  conditional  interval 
failure  probabilities  obtained  by  the  actuarial  method.) 
Moments  and  percentiles  of  F(t)  and  F^(t)  are  easily  obtained. 
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3.  SPECULATION  ABOUT  THE  PROPERTIES  OP 
THE  MAXIMUM  LIKELIHOOD  ESTIMATOR 

The  only  property  that  can  be  claimed  for  the  estimator 
F(t)  is  that  it  is  a function  of  the  sufficient  statistics. 

That  is  a consequence  of  the  maximum  likelihood  method  used 
to  obtain  the  estimator,  (Lindgren  (2]  p.225  ) . However,  many 
other  desirable  large  sample  properties  can  be  attributed 
to  maximum  likelihood  estimators . Asymptotic  properties  of 
umbiasedness,  efficiency,  consistency  and  normality  of  fj  and 

A 

Pj,  j=l,2, . . . ,k+l,  follow  from  verification  of  regularity  condi- 
tions on  the  second  derivatives  of  the  distribution  F(t)  , 

(Wilks  ij]  chapter  12  ) • This  reduces  to  verification  of  the 

same  regularity  conditions  on  F^Ct)  because  d~Pi) 

i=l 

satisfies  the  conditions.  The  user  of  the  model  (2)  may 
assume  the  conditions  are  satisfied  by  Fj^(t),  and  the  large 
sample  properties  are  immediate  consequences. 
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The  sample  proportions  of  survivors  , j*l,2,...,k  are 
so  intuitively  appealing  that  one  may  expect  desirable  small 
sample  properties  usually  attributed  to  sample  proportions — 
unbiasedness,  efficiency,  and  perhaps  others.  No  investiga- 
tion of  this  possibility  has  been  made. 

The  Glivenko-Cantelli  theorem  (Gnedenko  t4]  p.391)  states 
that  the  empirical  distribution  function  is  consistent.  The 
estimators  f^  appear  to  be  the  jumps  of  the  empirical  distri- 
bution function  conditioned  on  the  probability  of  not  being 
removed  at  previous  inspections.  However,  that  probability 

A, 

is  not  known  but  estimated.  A proof  of  consistency  of  Fj^(t) 
would  have  to  deal  with  that  estimate  and  then  app''y  the 
Glivenko-Cantelli  theorem. 

If  Fj(t)  is  assumed  to  be  continuous,  the  absolute 
difference  between  the  empirical  distribution  function  and 
Fj(t)  has  known  asymptotic  distribution  (Billingsley  [5Jp.l04). 
This  result  of  Kolmogorov  may  be  extended  to  the  maxi- 

A 

mum  likelihood  estimator  Fj^(t)  of  F(t)  if  the  conditional 

y\ 

form  of  the  jumps  f j , j =1 , 2 , . . . ,k+l  is  exploited.  There 
has  been  no  attempt  to  do  this. 

/s 

It  would  be  of  some  interest  to  compare  F(t)  with  the 
empirical  distribution  function  (which  is  an  unbiased  func- 
tion of  the  maximum  likelihood  estimator  of  a distribution 
function).  The  efficiency  of  F(t)  is  probably  higher  but 
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perhaps  not  high  enough  to  justify  the  slight  increase  in 
computation  required. 

In  comparison  with  the  conditional  interval  failure  proba- 
bilities  (called  failure  rates)  estimated  by  the  actuarial 
method,  the  corresponding  estimators  derived  from  F(t) 
certainly  have  better  large  sample  properties  since  they 
are  not  constrained  to  have  constant  failure  rate  in  each 
of  a fixed  set  of  intervals.  Where  data  is  more  complete, 

A 

the  intervals  between  jumps  of  Fj^(t)  are  small,  giving  more 
precision  during  time  periods  when  removals  are  more  likely. 

A, 

Where  data  are  scarce  because  few  removals  occurred,  F(t) 
is  still  an  improvement  over  the  actuarial  estimate  which  is 
a linear  extrapolation. 

The  feature  of  different  interval  si::es  indicates  a 
computational  problem  that  may  become  serious.  The  esti- 

A 

mate  F(t)  jumps  at  each  observed  usage  removal  and  inspec- 
tion time,  and  as  data  increases,  some  kind  of  fixed  interval 
grid  may  have  to  be  established  on  [0,  in  order  to  keep 

the  estimator  in  a reasonable  amount  of  computer  storage.  This 
need  creates  a new  statistical  problem  of  optimal  interval 
size  and  location. 

There  is  some  rate  of  convergence  comparison  of  interval 
estimators  for  restricted  families  of  distributions  in  chapter  5 
of  a book  by  Barlow,  Bartholomew,  Bremner,  and  Brunk[6].  Similar 
comparison  may  be  possible  for  the  empirical  distribution 

A A 

function,  F(t),  F(t)  constrained  to  optimal  intervals,  and  the 
estimator  given  by  the  actuarial  method. 
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Before  any  great  effort  is  spent  proving  these  specula- 
tions,  simulation  should  be  used  to  compare  F](t)  with  the 
empirical,  distribution  function  and  the  actuarial  method. 

A 

If  the  comparison  is  favorable  to  Fj(t) , it  could  be  used 
while  its  more  subtle  qualities  are  examined. 


SECTION  III 


AN  HIERARCHY  OF  ENGINE  REMOVAL  TIME  DISTRIBUTIONS 


1.  DESCRIPTION  OF  THE  HIERARCHY 

An  hierarchy  can  be  constructed  of  distribution  functions 
of  engine  flying  hours  till  removal  for  each  engine  type.  In 
Part  1,  feunllles  of  distribution  functions  will  be  defined 
with  an  hierarchical  relationship  that  lends  itself  to  sequential 
likelihood  ratio  tests  to  determine  which  family  represents  engine 
life  data  for  an  engine  type.  The  point  of  this  construction  is 
that  estimators  of  distribution  functions  for  the  more  restric- 
tive families  are  more  accurate  in  a statistical  sense.  Also, 
updating  the  official  failure  rate  can  be  done  statistically 
once  the  family  representing  the  engine  life  distributions  has 
been  identified.  These  latter  topics  are  discussed  in  Parts 
2 and  3.  Part  4 anticipates  some  criticism. 


Let  F(t)  denote  the  cumulative  distribution  function  of  the 
flying  hours  till  removal  of  a particular  type  of  engine.  Script 
7 denotes  a feunily  of  distribution  functions.  The  feunilies 
proposed  for  the  hierarchy  are  listed  in  increasing  order  of 
restrictiveness . 

7 2 ~ all  absolutely  continuous,  discrete,  and 

mixed  distributions  with  support  [0,  t„.„l  , t„.„  <» 

7 jj  ~ all  distributions  of  form 


[t] 


F(t)  « 1-  F, (t)  n 
^ i-1 


(1-p^)  0£t£t. 


max 


where  Fj^(t),  the  tail  c.d.f.  of  Fj^(t),  is  assumed  to 
be  absolutely  continuous  and  i is  the  index  of 
the  last  inspection  prior  to  t,  i^^j  * 0,1,2,... 
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(The  product  (l-p^)  is  assumed  to  be  1 for 

i=l  ^ 

i[tj*0.)  The  discontinuity  representing  removal 
of  engines  at  t is  contained  in  p.  for  l^ir^.  i 

Then  F,  (t)  has  as  its  support  [0,  t„_„]  . 

X InclX 

K J1J  ~ Various  restricted  families  of  distributions  such 
as  IFR  (increasing  failure  rate) , IFRA  (increasing 
failure  rate  on  the  average)  [7] > NBU  (new  better 
than  used) , NBUE  (new  better  them  used  in  expected 
life)  [8],  DMRL  (decreasing  mean  residual  life)  [8] 
and  all  the  familiar  distributions  of  nonnegative 
random  variables, 

7 - F(t)  - 1-e"^^  (1-p)  0<p<l,  0<t<-* 

The  last  two  families  contain  distributions  for  which  optimal 
replacement  policies  are  known  [8].  Other  families  may  be  of 
interest  such  as 
^ 

7 ^ - F(t)  = l-Fj^(t)  (1-p)  0<t<tjjj^j^,  0<p<l 

The  hierarchical  relation  among  these  families  is 

7 V Z>  7 IV 

7 III  ^ 7 V 3 7iv 

(The  inclusion  7^^^  37^  depends  on  which  restrictions  in 
7 are  chosen.)  Some  hierarchical  relations  among  the 
restricted  families  in  7 have  also  been  established. 

* These  families  could  be  truncated  at  if  desired. 


7 1 3 7 II  3 -I 
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2.  h 8EQUEHTIAL  LIKELIHOOD  RATIO  TEST  FOR  THE  APPROPRIATE  FAMILY 
Since  the  hierarchy  is  constructed  so  that  there  is 
containment  and  since  maximum  likelihood  estimators  of  F(t) 
exist  for  each  family,  a sequence  of  likelihood  ratio  tests 
may  be  used  to  determine  which  feunily  describes  data  from  an 
engine  type.  (The  maximum  likelihood  estimators  for  7 and 
are  derived  in  Section  II.)  Strong  preference  for  family 
7 j j begs  a test  to  determine,  for  each  engine  type,  whether 
F(t)e  (Reasons  for  this  preference  are  discussed  in 

Section  II.)  This  is  accomplished  by  the  likelihood  ratio  test 
of 

Hq;  F(t)  e vs. 

H^:  F(t)  e 7 j 

The  procedure  is  to  accept  the  hypothesis  if  the  likelihood 


ratio 


x^  ^ L(Tj^,...,T^iF(-)e7ii) 


max 
F(*)e 


max 


L(Tj^,...,T^1f(-)e7j) 


(3) 


F(.)e  7j 

is  sufficiently  large.  (The  likelihood  function  L (*|*)  is 
based  on  a sample  engine  removal  tiroes  given  the 

family  containing  the  underlying  distribution  function  F (t) . ) 
If  H^  is  accepted,  then  proceed  to  test 
Hq:  F(t)e7  vs. 

"a=  r(t)c7  „ 

in  the  same  manner.  Continue  to  test  for  7 jy  end  7y  until 
the  null  hypothesis  is  rejected  or  7 is  accepted. 
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i 

The  result  is  that  the  family  7^,...,  7^  may  be 

Identified  for  each  engine  type,  and,  as  the  test  proceeds 
from  7^  through  7^,  estimators  of  the  underlying  life 
distribution  function  beccxne  successively  better  in  a statis- 
tical sense.  There  are  other  potential  benefits  for  identifying 
a family.  The  f2unlly  7 jy  of  exponential/geometric  distri- 
butions allows  analytical  solution  of  the  replacement  problem 
described  and  simulated  in  Section  IV.  Some  subsets  of  family 
7 III  allow  bounds  on  the  solution  of  the  replacement  problem 
simulated  in  Section  IV,  [8]. 

Improvement  over  the  actuarial  method  is  possible  for  all 
f2unilies  7^,...  7^  because  the  estimators  for  each  family 
are  not  restricted  to  fixed  intervals  as  is  the  actuarial 
estimator.  E.g.,  for  7^  the  roost  common  estimator  is  the 
empirical  d.f.  with  jumps  of  equal  size  at  each  observed 

engine  removal  time.  However,  for  large  samples  the  empirical 
distribution  function  requires  excessive  data  storage,  and  some 
consideration  will  have  to  be  given  to  interval  estimators,  but 
not  necessarily  to  equal  size  intervals  as  in  the  actuarial 
method. 
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3.  UfMTING  THE  OFFICIAL  FAILURE  RATE 

The  Aerospace  Engine  Life  Committee  periodically  decides 
whether  to  update  the  "Offical  Failure  Rate"  curve  by  visual 
comparison  of  the  Official  Failure  Rate  curve  and  the  failure 
rate  computed  by  the  actuarial  method  from  current  data,  be- 
cause data  from  current  engine  operation  should  be  incorporated 
if  it  will  add  to  the  accuracy  of  estimates.  Also,  data 
from  earlier  operations  should  be  eliminated  if  it  appears 
that  changes  in  maintenance,  engine  modifications,  and 
usage  have  resulted  in  a change  in  the  failure  rate  of  an 
engine  tyne.  The  decision  of  the  Aerospace  Engine  Life 
Committee  can  be  assisted  by  a statistical  comparison  within 
a family  of  distributions  for  each  engine  type. 


The  statistical  procedure  requires  several  comparisons  of 
data  from  one  quarter  with  data  from  several  other  quarters. 
Some  definitions  are  necessary  before  description  of  the  pro- 
cedure. Let  F^(t),  F2(t),...,  denote  the  true  underlying 

distribution  functions  of  engine  flying  hours  at  removal  in 
quarters  1,  2,...,k.  Quarter  1 is  the  first  quarter  for  which 
data  is  used  in  calculation  of  the  Official  Failure  Rate,  and 
quarter  k-1  is  the  last  quarter  of  data  which  is  used  in  the 
calculation.  Quarter  k is  the  current  quarter  with  engine 
removal  time  data  under  consideration  for  incorporation  into 
the  Official  Failure  Rate.  Let  the  underlying 

c.  d.  f.  of  engine  removal  times  throughout  quarters 

j,  j+1, . . . ,k-l;  j=l,2, . . . ,k-l.  The  latter  may  not  be  the  same 
as  Fj^  (t) , . . . (t)  and  Fj^(t)  due  to  changes  in  maintenance, 

usage,  or  design  of  the  engines. 
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At  the  time  of  data  analysis  each  quarter,  a set  of 
hypotheses  should  be  tested  as  shown  in  Figure  2.  An 
explanation  of  the  procedure  may  be  helpful.  First  test 
whether  to  incorporate  the  new  data. 

“o*  “ iVl^^^ 

Fjj(t)  ^ te[0, 

If  is  accepted,  then  data  from  the  current  quarter  should 
be  considered  for  incorpor^tion  into  the  Official  Failure  Rate 
but  not  before  this  next  test  to  determine  whether  the  oldest 
data  continue  to  be  used. 

H3:  Fj^(t)  te[0, 

If  hypothesis  H2  is  accepted,  the  estimator  of  engine  removal 
time  distribution  should  be  revised  to  include  data  from  all 
k quarters. 

If  either  hypothesis  rejected,  additional  study  is 
necessary  to  determine  which  data  should  be  used  to  estimate 
the  engine  removal  time  distribution  and  Official  Failure 
Rate.  If  is  rejected,  the  data  for  quarter  k seems  to  have 
come  from  a population  of  engines  that  does  not  resemble 
engines  in  quarters  l,2,...,k-l  in  its  flying  hour  1 <:e. 
Further  tests  should  be  made  to  compare  F.  (t)  with  (t) 


If  passes  statistical  comparison  with  some 

j,  the  official  engine  life  c.d.f.  should  be  updated  to  an 

A 

estimate  jFj^  including  data  from  the  current  quarter.  If 
not,  then  the  official  engine  life  data  should  be  based  on 
the  current  quarter,  but  not  until  there  is  an  explanation 
for  the  apparent  change  and  an  assessment  of  whether  the 
current  quarter  usage  is  representative  of  the  future.  This 
decision  should  remain  the  responsibility  of  the  Aerospace 
Engine  Life  Committee.  It  should  also  be  recognized  that 
all  other  statistical  conclusions  in  this  section  are  only 
recommendations.  These  recommendations  will  be  accompanied 
by  a statistical  level  of  confidence,  a measure  of  the 
possible  error  in  the  recommendation. 

As  indicated  in  Part  2 of  this  section,  increased 
statistical  accuracy  of  the  recommendations  is  the  result  of 
the  hierarchical  relationship  among  families  T j,...,  Ty 
Once  the  family  has  been  established  as  in  Part  2,  comparison 
of  distribution  functions  within  each  family  as  required  in 
this  part  may  be  carried  out  with  the  aid  of  existing  statistical 
tests . For  example : 

^ j.:  Kolmogorov  - Smirnov  test, 

"y  Kolmogorov  - Smirnov  test  on  Fj^(*)  and 

multinominal  test  on 


111’  Likelihood  *ratio  tests  which  exist  for 

certain  restricted  families  (Barlow  [7] ) 
Tests  of  the  exponential  failure  rate  r 
the  survival  probability  p,  and 
y Same  as 
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4.  CAUTION  ON  THE  USES  OF  THIS  SECTION 

The  major  limitation  is  that  further  development  is  needed. 
This  development  should  proceed  in  conjunction  with  the  use  of 
the  procedures  in  this  section  in  management  decision  making. 
Statistics  can  only  specify  the  probability  of  error  in  an 
hypothesis  test.  Management  decision  makers  needs  to  deter- 
mine what  tests  are  necessary  and  c.n  acceptable  levtl  of  risk. 
This  section  is  proposed  in  anticipation  of  decision  making 
needs  and  for  use  in  obtaining  better  estimates  within  the 
existing  data  structure. 

Fortunately,  most  of  the  statistical  work  remaining 
to  be  done  is  development,  rather  than  high  risk  basic  research. 
Maximum  likelihood  estimators  of  F(t)  exist  for  all  families. 
(See  Section  II  for  estimators  of  n v' 

restricted  distributions  in  jy  also  have  estimators.)  The 
remaining  work  is  to  develop  sequential  likelihood  ratio  tests 
and  other  comparison  tests  that  are  convenient  for  use. 
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SECTION  IV 


FORMULATION  AND  SIMULATION  OF  THE 
ENGINE  REPLACEMENT  PROBLEM 


1 . INTRODUCTION 

Stripped  of  all  complications,  engine  inventory  manage- 
ment requires  reorder  quantities  of  each  engine  type  each 
quarter  based  on  the  flying  hour  program  for  that  quarter. 

The  reorder  quantity  is  equal  to  the  number  of  replacements 
taken  from  inventory.  The  number  of  replacements  is  the  number 
of  engines  that  are  removed  but  not  reinstalled,  so  a predic- 
tion of  the  number  of  removals  as  a function  of  the  flying 
program  is  required  for  inventory  management.  That  is  the 
engin*  replacement  problem  to  be  studied  in  this  section. 

Even  under  the  simplification  that  there  is  only  one 
single-engine  plane  at  a base,  the  distribution  of  the  number 
of  removals  cannot  be  simply  computed  as  a function  of  the 
removal  time  distribution  and  the  flying  program.  The  reason 
is  that  the  replacement  engines  are  not  necessarily  new  but 
have  already  accumulated  a known  amount  of  flying  time.  The 
only  exception  occurs  when  the  removal  time  distribution  has  the 
"memory less  form  of  the  exponential  and  geometric  distributions, 

F(t)  = l-e“^^(l-p)  X,  t>0,  0^<1  (4) 

where  i^^j denotes  the  index  of  the  last  inspection,  if  any, 

prior  to  or  at  time  t.  In  this  case  the  replacement  process  becomes 
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a renewal  process  for  which  many  results  are  known  (Feller 
[ 9]  and  Cox  [ 10]  ) . 

Despite  obstacles  to  solution,  formulation  of  the  single 
engine  replacement  problem  is  of  value  to  point  out  analyti- 
cal limitations,  to  suggest  approximations,  and  to  identify  the 
mathematical  model  that  can  be  simulated.  The  single  engine 
simulation  can  be  adapted  to  multi-plane,  multi-engine  simula- 
tion with  complications  sufficient  to  defy  analysis. 

Computer  models  for  simulation  of  the  engine  replacement 
problem  are  given  in  Part  3 for  single  quarter  and  multiple 
quarter  planning  periods.  If  the  sequence  of  replacement 
engines  (with  known  accumulated  flying  hours)  is  random,  the 
conditional  distribution  of  the  nvunber  of  engines  in  stock 
with  less  than  a certain  age, given  the  number  of  replacements , 
has  binomial  form.  This  distribution  may  be  used  to  simulate 
a random  replacement  sequence.  Further  combinatorial  analysis 
may  yield  this  distribution  if  there  are  many  operating 
engines.  However,  engine  replacement  may  be  simulated 
regardless  of  pending  analysis,  and  the  event  by  event  type 
computer  model  proposed  is  much  quicker  in  execution  than 
interval  by  interval  simulation  based  on  the  actuarial 
method. 

Another  reason  for  the  event  type  computer  model  is  that 
it  may  be  easily  modified  to  incorporate  antithetic  random 
variables.  These  negatively  correlated  random  variables  are 
used  to  reduce  the  variance  of  means  and  percentiles  in  simu- 
lation of  waiting  times  in  single  server  queuing  systems 
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(Mitchell  [ll3  ) . The  same  technique  and  proof  can  be  applied  to 
simulation  of  the  mean  value  function  of  a renewal  process. 
Although  data  are  not  very  encouraging,  the  use  of  antithetic 
variates  may  reduce  the  variance  of  the  sample  mean  in  simula- 
tion of  the  engine  replacement  problem. 

Further  analytical  and  empirical  study  of  the  engine 
replacement  prediction  problem  is  recommended  in  conjunction 
with  analysis  of  the  buffer  quantity  to  be  reordered  each 
quarter.  Approximations  should  be  compared  since  some  are 
being  used  without  any  current  justification.  The  event  by 
event  simulation  technique  should  be  incorporated  since  it 
saves  computer  time  in  what  promises  to  be  a large  scale 
simulation.  Also  there  should  be  further  investigation  of 
the  application  of  antithetic  variates  to  simulation  of  the 
number  of  engine  replacements. 

2.  FORMULATION  OF  THE  ENGINE  REPLACEMENT  PROBLEM 

Fundamental  to  the  engine  replacement  problem  is  the 
random  process  N(t)  representing  the  number  of  engine  replace- 
ments on  a single  engine  plane  necessary  to  complete  t flying 
hours.  (Modifications  can  be  made  for  a multi-engine  problem 
with  a random  sequence  of  replacements.)  The  distribution 
F(t),  t^O,  of  flying  hours  between  removals  is  assumed  to  be 
known.  The  currently  installed  engine  has  already 
accumulated  t^  flying  hours  and  the  replacement  engines  have 
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•<•*  flying  hours.  They  replace  the  original 
engine  in  that  order.  Denote  their  residual  lives  given 
ages  tg,  t^,...tg  by  R(tQ) , . . . ,R(t^) . These  random  variables 
are  assumed  independent  and  have  the  distribution 

F(t^+v)-F(t^) 


PtR(tj^)  < v|removal  time  > tj|^] 
for  v>0,  and  i«0,l,2, . . . ,s. 


l-F(t^) 


(5) 


The  relation  between  partial  sums  of  residual  lives  and 
the  number  of  replacements  in  a time  interval  t may  be  used 
to  obtain  the  distribution  of  N (t) . If  an  engine  is  not 
removed  in  10, tl,  then  R(tQ)>t,  and  no  replacement  occurs.  If  HCt^) 
et0,t]  and  the  second  engine  survives,  then  R(tQ)+R(tj^)>t,  and 
only  one  replacement  occurs.  The  general  relation  between 
N(t)  and  partial  stuns  of  R(tj^)  may  be  expressed  in  terms  of 
the  equivalent  events: 


RCtg) >t 

R(tQ) <t<R(tQ)+R(t^) 


(6) 


N(t)  = n 


n-1  n 

I R(t.)<t<  I R(t.) 
^ - i=o  ^ 


Therefore,  the  probability  distribution  of  N(t)  is  known  if 
n-1  n 

Ft  I R(t.)<t<  I R(t.)J  = P[N(t)=nJ 
i=0  ^ i=0  ^ 

can  be  determined  for  all  n.  That  requires  the  distribution 
of  partial  sums  of  residual  lives  given  the  ages. 
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J R ( t • ) / n*l f 2 f f 

i-0  ^ 


The  density  of  the  sum  of  two  Independent  non- 
negative random  variables  is  obtained  by  the  convolution 
formula 


^PtXi+Xj<tl 


(t-u)dP2  (u) 


(7) 


where  is  the  density  of  and  F2(*)  is  the 

c.  d.  f.  of  X2*  If  X^  and  X2  have  the  same  distribution,  the 
convolution  is  of  simpler  form,  but  it  is  still  inconvenient 
to  compute  except  for  special  cases.  If  X^  and  X2  are  two 
residual  lives,  Equation  (7)  becomes  the  convolution  of  two 
distributions  (5)  with  possibly  different  ages.  There  is  no 
computationally  convenient  representation  for  such  convolutions 
even  though  the  integral  formula  is  easily  expressed  as 


^ PlR(t^)  + R(tj^)  <t] 


f (tjj+U) 
l-FCtj)^  . 


f (t-U+tj^) 
l-F(tj^) 


du 


where  f(«)  denotes  the  density  of  ?{•),  assumed  to  exist  for 
this  illustration. 


Several  approximations  may  be  adequate  for  determination 
of  the  reorder  quantity  in  certain  situations.  One  of  them 
is  in  current  use,  and  the  others  warrant  investigation  also. 

1)  Assume  that  every  engine  lasts  its  expected 

residual  life,  and  compute  the  corresponding  N(t). 
If  the  engine  ages  and  the  sequence  of  replacements 
are  known,  then  N(t)  is  determined.  This  method  is 
now  used  with  some  modifications. 
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2)  If  the  flying  program  represented  by  t is  small 
relative  to  the  residual  life  expected  for  the 
currently  installed  engine,  an  adequate  approxi- 
mation to  the  distribution  of  N(t)  may  be  obtained 
by  numerical  convolution  for  the  first  few  values 
of  P[N (t)*n] . 

3)  If  the  value  of  t is  large  relative  to  the  residual 
lives  of  the  first  engine  and  its  replacements, 

the  central  limit  theorem  for  approximately 
identically  distributed  random  variables  (Gnedenko  [ 12  ] 
chapter  8 ) niay  yield  an  adequate  normal  approxi- 
mation for  the  partial  sxims  of  residual  lives. 

4)  If  the  distribution  of  removal  times  F(t)  is 
sufficiently  similar  to  the  exponential  and 
geometric  distributions  (4 ) , it  may  be  practically 
memoryless.  The  replacement  process  then  becomes 
a renewal  process.  Analysis  of  exponential  and 
geometric  renewal  processes  is  essentially 
complete . 

5)  If  instead  of  using  the  known  accximulated  flying 
hours  for  each  replacement  engine,  the  equilibrium 
age  distribution  is  assumed  for  each  replacement, 
then  the  generating  function  of  the  Laplace  transform 
of  the  distribution  of  N(t) 

S 00 

G(z,H>)  = I f e~'^^  P[N(t)»n]dt 


for  ¥^0,  0<z<l  may  be  found.  The  formula  for  G(z,*f)  is 
given  by  Cox  [10]  on  page  38.  If  this  yields  an  adequate 
approximation  to  the  distribution  of  N(t),  then  it  may 
not  be  necessary  to  incorporate  the  actual  ages  into 
prediction  of  N(t).  However  there  is  a notorious  prob- 
lem of  inverting  the  generating  function  of  a Laplace 
transform  to  obtain  the  distribution  function  of  a random 
process. 

6)  Simulation  of  values  of  N(t)  can  be  done  given  the 

distribution  F(t)  of  removal  times,  ages  t^,  t^,...,t^, 
the  residual  life  distribution  ( 5) , and  the  relations 
(3) . This  method  is  strongly  recommended.  It  will 
be  described  in  part  3 of  this  section. 

Given  the  distributions  of  N(t)  conditioned  on  all  possible 
sequences  of  replacements  t. ,...,t  , the  distribution  of  N(t) 

X 9 

for  a random  replacement  policy  may  be  determined.  An 
ingredient  necessary  for  calculating  the  probability  of  every 
possible  replacement  sequence  is  the  age  distribution  of  the 
inventory  after  a specified  nusiber  of  replacements  have 
occurred.  Each  replacement  is  a random  selection  from  the 
inventory  with  age  distribution  which  is  given  next. 

Let  DQ(x;t)  be  a step  function  that  tells  the  proportion 
of  engines  initially  in  stock  with  accumulated  flying  hours 
less  than  or  equal  to  t,  x»0,l,2, . . . ,s,  t^O.  D^(x;t)  gives 
the  same  information  after  the  first  replacement  has  been  made, 
except  that  x>x(t)  is  the  value  of  a random  function  of  t 
because  it  is  not  known  which  engine  is  chosen  from  stock. 

4l 


D^(x;t)  is  the  probability  function  after  the  yth  replacement 
y*0»l, . . . ,8<». 


Since  the  original  inventory  of  engines  had  accumulated  flying 

hours  t^,  t2r.».»tgr  the  probability  function  D^(x;t)  will  change 

only  at  those  times.  Let  t represent  the  time  value  t^e[t. ; 

i»l,2,...,s]  for  which  exactly  x engines  in  the  original  stock 

had  t or  fewer  accumulated  flying  hours.  Then  D (x,t)  can  be 
X y 

represented  in  binomial  form 


Dy(x-z;  t^)  . (J)  (|)*  (1-1)^=^ 

for  i>l,2,...,s}  and  z>0,l,2, . . . ,min(x,y)  . A few 

examples  may  be  sufficient  to  convince  the  skeptic.  Let 
tj^£t2ji. . .^tg  and  examine  the  function  Dj^(x;t)  for  the  engine 
age  distribution  after  the  first  replacement  is  installed: 


D^(0;t^)  s ^ is  the  probability  the  newest  engine  was 
installed, 

Di(l;ti)  * 1-^  if  the  newest  engine  was  not  installed; 
all  other  Dj^(x;  tj^)  are  zero.  The  general  form  of  Dj^(x;t)  is: 


D, (x-1;  t ) = x/s  if  an  engine  of  age  <t  was  chosen  as 
replacement,  and 

Di(x;  t^)  = 1-x/s  if  not. 

Additional  combinatorial  work  will  be  necessary  to  obtain 
the  distribution  of  the  number  of  replacements  for  a fleet 
with  more  than  one  engine  in  operational  status  supplied 
from  the  same  inventory  of  spares.  The  formidable  nature  of  the 
problem  is  sufficient  reason  to  suggest  simulation. 

42 


iVi 


3. 


SIMULATION  OF  THE  ENGINE  REPLACEMENT  SYSTEM 


It  is  necessary  to  simulate  the  engine  replacement  sys- 
tem in  order  to  obtain  some  information  about  the  distribu- 
tion of  the  number  of  replacements  for  a specified  flying  pro- 
gram. In  the  absence  of  any  solution  to  the  problem  formulated 
in  part  2 of  this  section,  simulation  will  provide  data  for  com- 
parison with  approximations  and  bounds,  and  it  will  also  provide 
a prediction  method  that  can  be  incorporated  into  engine  manage- 
ment procedures. 

The  first  flow  chart.  Figure  3,  outlines  a program  to  simu- 
late operation  of  a single-engine  plane  for  a specified  number 
of  flying  hours  and  a specified  sequence  of  replacement  engines 
with  known  accumulated  operating  times.  Incorporating  random 
selection  from  the  stock  of  replacement  engines  is  easily  done. 
Simulation  of  a multi-engine  fleet  is  next  illustrated.  It  is 
assumed  that  each  operating  engine  must  fly  an  equal  part  of 
the  flying  program,  but  this  assumption  is  easily  removed.  The 
flying  hour  program  for  each  engine  type  is  given  on  a quarterly 
basis.  For  long  range  planning,  a simulation  is  necessary  for 
many  quarters.  The  second  flow  chart.  Figure  4,  describes  an 
event  by  event  type  of  simulation  that  is  being  incorporated  into 
a master  program  by  L-  A.  Coco  of  AFLC/MMOMA. 

The  program  in  Figure  1 requires  as  input  accumulated 
flying  times  for  the  installed  engine  t^  and  for  all  spares 
ti,...,ts»  the  removal  time  c.  d.  f.  F(t),  the 
flight  plan  D,  and  the  run  size  K.  Each  cycle  simulates  the 


flying  program  and  generates  a value,  N(D),  for  the  number  of 
replacements  required  to  achieve  the  flying  progreun.  Resi- 
dual lives  R(t^)  are  generated  from  the  distribution (5)  by 
any  method,  (see  Fishman  [13]  chapter  8).  Statistics  computed 
on  the  output  values  of  N(0)  depend  on  the  application,  but 
there  are  standard  computer  programs  for  most  needs.  If  de- 
sired, flight  plan  D and  ages  tg,  tj^,...,tg  could  be  randomly 
generated  as  part  of  the  program.  A random  choice  of  replace- 
ment engines  from  stock  is  easily  incorporated. 

A more  complicated  program  is  required  to  simulate  multi- 
period operation  of  a multi-engine  fleet.  Flexibility  is  provided 
to  allow  incorporation  of  a different  flight  plan  and  changing 
inventory  and  age  distribution  for  each  period.  Each  cycle 
through  the  program  generates  the  number  of  replacements  re- 
quired to  satisfy  the  flight  plan  for  each  period.  The  flow- 
chart in  Figure  4 shows  just  one  such  cycle.  The  master  program 
which  will  use  Figure  4 generates  transfers  to  and  from  the 
inventory  and  ages  of  engines  put  into  inventory. 

Variables  used  in  the  program  are  defined  as  follows: 

N = the  number  of  quarters  to  be  simulated; 

PROG(I)  = flying  program  for  quarter,  1=1, 2 , . . . ,N; 

AGE  (IN)  = age  of  installed  engine,  IN  = 1,2,...,  total 
numbers  of  engines  installed; 

SAGE(JN)  = age  of  spare  engine,  JN  = 1,2,...,  upper  bound 
on  the  number  of  spares  available; 

RES (IN)  5 residual  life  of  an  installed  engine  after 
completion  of  quarterly  flying  program; 
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EPROG(I)  = flying  hour  program  for  each  engine  and  quarter; 

RES  H residual  life  of  an  engine  (RES  t AGE (IN)  tells 
the  total  number  of  flying  hours  on  an  engine 
before  removal) ; 

INSTIN(I)  i the  number  of  installed  engines  in  quarter  1; 

INV(I)  H the  number  of  spares  in  quarter  I. 

The  routine  in  the  dotted  box  is  provided  by  the  master  program. 
A varying  maximum  operating  time  for  engines  may  be  incorporated 
into  the  program. 

The  event  by  event  type  of  simulation  will  sav'’!  execution 
time  compared  with  the  interval  by  interval  simulation  of  each 
engine  now  used  in  the  master  program.  In  that  program,  each 
quarter  is  divided  into  many  equal  size  intervals  and  each 
engine  is  simulated  to  see  whether  it  survives  successive  in- 
tervals. While  this  method  conveniently  incorporates  actuarial 
failure  rates,  it  is  not  necessary.  Survival  probabilities 
obtained  by  the  actuarial  method  may  be  converted  into  the 
c.  d.  f.  F(t)  required  as  input  in  Figure  4.  An 

alternative  is  to  use  the  estimator  derived  in  Section  V. 

4.  APPLICATION  OF  ANTITHETIC  VARIATES  TO 

THE  REPLACEMENT  PROBLEM  SIMULATION 

The  complicated  program  required  to  simulate  engine  manage- 
ment systems  may  strain  the  capacity  of  available  computers. 

One  method  to  reduce  variability  of  output  for  a fixed  run 
size  or  to  reduce  run  size  for  fixed  level  of  variability  is 
the  method  of  antithetic  variates  (Hammersley  and  Handscomb 
[14]).  This  method  has  long  been  known  in  statistical  sample 
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design  but  has  received  little  use  in  simulation.  The  method 
will  be  illustrated  in  simulation  of  the  mean  value  of  ran- 
dom variable.  Application  of  the  technique  to  simulation  of  the 
replacement  problem  of  part  3 is  easy  to  do.  It  has  been  tested 
for  a renewal  and  a replacement  process  with  uniform  and  expon- 
ential removal  times  with  the  same  mean.  In  most  tests, 
anitthetic  variates  gave  a reduction  in  the  variance  of  the 
mean  number  of  removals.  An  outline  of  a proof  is  given  that  this 
is  true  for  all  distributions  when  simulating  the  renewal  process 
and  the  replacement  process. 

The  simple  program  flowchart  in  Figure  5 is  designed  to 
estimate  the  expected  value  of  a random  variable  with 

_ N 

c.d.f.  F(x).  The  sample  mean  X = Z Xi/N  is  an  un- 

i=l 

biased,  consistent,  and  otherwise  well  behaved  estimator  of 
the  expected  value  of  a random  variable.  Its  variance  is 

_ ( N ) 2 

Var  X = < E Var  X^  + 2 Z Cov  (Xi,X^)  > /N  (8) 

(i=l  i<j  ^ ) 

The  covariance  term  enters  into  the  calculation  if  the  random 
variables  X^^,...,Xjj  are  dependent.  The  "trick"  of  the 
antithetic  variate  method  is  to  generate  the  sequence  so  that 
Cov  (Xj^,Xj)  ^ 0 for  all  i and  j.  In  some  cases,  exponential 
and  uniform,  the  use  of  negatively  correlated  random  number 
sequences  for  generation  of  Xi,...X„  will  yield  negatively 
correlated  X^  (Fishman  [13]  p.  320)  . 
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INPUT  RUN  SIZE  N 
AND  DISTRIBUTION  FUNCTION  F(x) 


J*1 ,s*o 


GENERATE  RANDOM 
NUMBER  U(J) 


GENERATE  RANDOM  VARI- 
ABLE X(J)  FROM  U(J) 


Figure  5.  Simulation  of  the  Expected 
Value  of  a Random  Variable 


49 


Two  methods  for  creating  negative  correlation  among  random 
numbers  are  to  use  the  complement,  1-random  number,  and  to 
generate  half  the  required  set  of  random  numbers  and  reverse 
the  sequence  for  use  in  the  second  half.  The  method  of  com- 
plements flowchart  is  given  in  Figure  6.  Both  methods  were 
tested,  and  the  former  seems  more  promising  for  simulation 
of  renewal  and  replacement  processes.  One  reason  for  this 
preference  is  that  many  random  numbers  must  be  generated  to 
give  one  value  of  the  number  of  removals  according  to  the 
relations  (6) . Reversing  the  sequence  of  random  numbers  used 
to  generate  the  R(tj^)  in  the  next  run  may  give  independent 
values  of  N(t)  if  n is  small  relative  to  the  number  of  random 
numbers  generated.  For  example,  suppose  a test  program 
generates  50  random  numbers  for  each  run,  U(l),  U (2) , . . .U (50) , 
which  are  used  to  generate  50  residual  lives.  If  t is  small, 
then  n may  also  be  small,  say  5.  On  the  second  run,  the 
sequence  of  random  numbers  will  bo  used  in  reverse  order, 

U(50),  U (49) , . . . ,U(1) , and  N(t)  will  probably  involve  only 
the  first  few  values  of  the  reversed  sequence,  independent 
of  U(l)  , U(2)  , . . , ,U(5)  . 


Figure  6.  Simulation  of  a Uniformly  Distributed 
Renewal  Process  with  Antithetic  Variates 


The  proof  of  the  variance  reduction  property  of 
complementary  antithetic  variates  will  only  be  sketched 
because  a search  is  now  under  way  for  a more  elegant  proof. 

The  proof  outlined  here  applies  the  results  of  Mitchell,  [11]. 
The  statement  to  be  proved  is: 

The  use  of  complimentary  antithetic  variates  reduces  the 
variance  of  the  sample  mean  in  the  simulation  of  the  mean 
value  function  E N(t)  of  a renewal  process  or  a replace 
ment  process. 

The  method  of  proof  is  to  use  Michell's  propositions  P1-P5 
to  show  negative  correlation  between  counting  functions  in 
antithetic  simulations  of  a replacement  process  for  a fixed 
length  of  time  t. 

The  method  of  proof  is  to  construct  a sequence  of  random 
variables  {Nj(t);  j^l}  for  fixed  t that  satisfy  Mitchell's 
propositions.  The  propositions  are  slightly  restated  here 
for  brevity  and  for  this  application. 

PI.  Nj^(t)  = 0 with  probability  one. 

P2.  There  exists  a sequence  of  functions  { f„}  such 
that  f^  maps  the  unit  n dimensional  hypercube 
into  the  non-negative  real  line. 

P3.  Nj  (t)  is  independent  of  Uj^  for  i> j . 

P4.  f^  is  non-increasing  in  for  all  j, 
l^j  ^n . 

P5.  The  stationary  distribution  of  N.  exists. 


Let  be  a sequence  of  independent,  identically 

distributed  non-negative  random  variables  generated  from  a 

c.d.f.  F(x)  by  inverse  transformations  of  a sequence  of 

uniform  random  variables  (U.)  . Let  (t)  be  the  function 

1 D 

that  counts  renewals  in  the  following  way 


Ni(t)  = 

0 

Njit)  =! 

0 

if 

^1 

> 

t 

1 

1 

if 

^1 

< 

t 

0 

if 

> 

t 

N3(t) 

1 

1 

if 

Xl 

< 

t n Xj^+X2>t 

2 

if 

^1 

+ 

X2  < t 

0 

if 

^1 

> 

t 

N (t)  = 

1 

if 

^1 

< 

t n x^^  + X2  > t 

4 J 

2 

if 

^1 

+ 

X2  < tnxi  + X2  + X3 

'3 

if 

^1 

+ 

X-  + X,  < t 

2 3 - 

etc . 


By  construction,  the  sequence{N^ (t)}  satisfies  PI, 

P2,  and  P3.  Since  the  Xj  are  generated  by  the  inverse 
transformation  F“^(Uj),  non-decreasing  in  , the  Nj (t) 
are  non-increasing  in  U^,  i ^ j.  In  the  context  of  the 
engine  replacement  process,  longer  engine  lives  (larger 
X^)  require  fewer  replacements  (Nj (t)  . Proposition  P5, 
the  existence  of  a stationarity  distribution,  is  satisfied 
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because  (t)  converges  to  the  renewal  or  replacement  counting 
function  N(t)  as  j increases.  The  (t)  represent  the  number 
of  replacements  if  there  are  only  j spare  engines,  and 
N(t)  is  the  number  of  replacements  if  there  are  unlimited 
spares.  Satisfaction  of  propositions  PI  - P5  is  sufficient 
for  proof  of  the  variance  reduction  property. 


SECTION  V 


AN  ANNOTATED  BIBLIOGRAPHY  OF  SOME 
POTENTIALLY  USEFUL  REFERENCES 

1.  CRITERIA  FOR  INCLUSION  IN  THE  BIBLIOGRAPHY. 

The  Articles  and  Research  Reports  described  here  were 
selected  for  their  relevance  to  three  levels  of  engine 
nan..gc  ^lent. 

1) i  Estimation  of  the  distribution  function  of  flying 

hours  to  engine  replacement  and  the  related  prob- 
lem of  predicting  the  number  of  replacements  re- 
quired to  meet  the  flying  program.  (These  are  the 
problems  studied  in  Sections  II,  III,  and  IV.) 

2)  The  logistics  problems  of  inventory,  repair,  dis- 
tribution, and  replacement  of  engines  of  a parti- 
cular type. 

3)  The  Air  Force  wide  engine  management  policies  and 
planning . 

General  results  and  applications  c re  balanced  to  give  a 
perspective  for  the  future.  More  recent  results  are  cited 
unless  there  appears  to  be  potentially  useful  work  that  has 
not  been  applied  to  engine  management. 

This  bibliography  is  not  claimed  to  be  comprehensive. 

It  does  not  include  applications  of  traditional  statistics 
to  stochastic  processes,  and  several  potentially  useful 
statistical  techniques,  such  as  time  series  analysis  and  spectral 
analysis,  have  been  omitted.  This  is  because  the  author 
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feels  that  their  potential  use  lies  furthc  the  future 

than  the  use  of  the  articles  cited  here  and  the  results  of 
Sections  II,  III,  and  IV.  Also  deliberately  omitted  is  the 
literature  on  system  reliability  under  assumptions  about  the 
structure  of  the  system  or  the  underlying  distribution 
function  of  lifetimes.  These  results  are  fairly  well  recognized, 
so  there  is  little  need  for  citation.  On  the  other  hand, 
this  Bibliography  includes  articles  from  1974  Department  of 
Defense  bibliographic  searches  on  maintainability  and  applica- 
tions of  stochastic  processes  to  engine  management,  so  it 
is  a fairly  thorough  indication  of  the  state  of  the  art 
familiar  to  the  armed  forces  and  deemed  relevant  by  the 
author. 
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2 . GENERAL  REFERENCES 

In  this  part  are  collected  general  references  of  potential 
use  in  engine  management  problems.  Two  specific  problems  were 
in  mind  during  the  selection  of  these  articles,  estimation  of 
engine  flying  hours  till  removal  and  prediction  of  the  number 
of  replacements  required  to  meet  the  flying  program.  However, 
the  generality  of  the  references  means  that  other  applications 
are  possible,  if  not  likely.  First  cited  are  articles  on  the 
estimation  of  stochastic  processes  under  some  assumptions  about 
the  processes.  Second  are  articles  that  may  collectively 
be  described  as  curve  fitting  techniques  with  potential  for 
applications  to  stochastic  processes. 

Consider  one  engine  mounting  position  on  an  operational 
aircraft.  A sequence  of  engine  replacements  occurs  in  that 
position,  and  data  for  each  engine  is  recorded  about  the  flying 
hours  till  removal  and  the  removal  cause.  If  all  replacement 
engines  are  new,  or  as  good  as  new,  the  stochastic  process 
counting  replacements  as  a function  of  flying  hours  may  be 
a "renewal  process".  If  replacement  engines  have  already 
accumulated  some  flying  time,  the  stochastic  process  will 
be  referred  to  as  a replacement  process.  There  are  several 
techniques  for  estimating  characteristics  of  renewal  processes 
that  are  potentially  useful  for  replacement  processes.  The 
references  focus  on  estimation  from  data  on  the  number  of 
replacements.  This  is  in  contrast  to  the  actuarial  technique 
and  the  estimator  derived  in  Section  II  which  estimate 


flying  hours  till  replacement  directly  from  flying  hours 
data.  It  is  possible  to  estimate  the  flying  hours  distribution 
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from  replacement  counts,  and  coljLection  of  replacement 
counts  is  easier  than  collection  of  flying  hours  data. 

A major  figure  in  the  development  of  statistics  for  re- 
newal processes  is  D.  R.  Cox.  His  article  ”0n  the  Estimation 
of  the  Intensity  Function  of  a Stationary  Point  Process", 

1965  [15],  deserves  citation  because  of  its  potential  for 
application.  The  article  provides  an  estimator  of  the  renewal 
intensity  function  which  can  be  integrated  to  obtain  an  esti- 
mate of  the  expected  number  of  replacements  for  any  number 
of  flying  hours. 

The  article,  "The  Determination  of  the  Distribution  Func- 
tion of  a Renewal  Process  from  a Single  Realization  Observed 
with  Some  Error"  by  P.  B^rtfai,  1966  [16],  is  cited  because  of  the 
interesting  possibility  of  providing  estimates  from  data  that 
contains  errors. 

Two  reference  books  follow  the  early  article  by  Cox. 

First  published  was  The  Statistical  Analysis  of  Series  of 
Events  by  D.  R.  Cox  and  P.  A,  W.  Lewis,  1966  [171.  Later, 

Lewis  was  editor  of  articles  presented  at  a symposium.  The 
title  of  that  collection  is  Stochastic  Point  Processes, 

1972  [18].  It  contains  some  articles  of  statistical  nature 
that  may  be  useful  in  engine  management.  In  "Multivariate 
Point  Processes",  1972  [19],  Cox  and  Lewis  provide  means  for 
characterizing  certain  kinds  of  multivariate  point  processes 
and  for  estimating  them.  This  last  article  may  be  useful  if 
it  is  necessary  to  analyze  tne  engine  replacement  process  as 
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a multivariate  process.  For  Instance  such  questions  as 
the  relation  between  calendar  age  and  flying  hours  and 
the  relation  between  flying  hours  and  maintenance  level 
require  bivariate  techniques. 

Until  this  point  in  the  bibliography,  methods  of 
estimation  that  make  assumptions  about  the  underlying  dis- 
tribution functions  have  been  avoided.  It  has  been  observed 
that  the  effect  of  periodic  inspection  and  other  factors 
result  in  an  engine  flying  hour  distribution  of  sufficiently 
complicated  nature  to  be  excluded  from  commonly  estimated 
faunilies  of  distributions.  However,  there  is  another  ap- 
proach which  may  be  of  some  use,  and  not  without  precedent 
in  estimation  of  engine  removal  rates.  That  approach  is  to 
make  some  assumptions  about  the  smoothness  of  the  distribu- 
tion function  or  its  associated  failure  rate  function.  (One 
precedent  is  that  the  engine  failure  rate  data  now  calcu- 
lated is  smoothed  by  moving  average,  extrapolated  by  linear 
regression,  and  then  presented  to  the  Aerospace  Engine  Life 
Committee  in  graphical  form  with  data  points  connected  by 
hand.)  There  is  no  functional  relation  between  the  smooth- 
ness assumptions  and  the  underlying  replacement  process.  It 
is  a matter  of  satisfaction  and  convenience  for  calculation. 

S.  A.  Krane  presented  the  article  "Analysis  of  Survival 
Data  by  Regression  Techniques",  1963  [20] , which  described 
estimation  of  the  integral  of  the  failure  rate  function  as 
a polynomial.  The  technique  is  for  large  seunples  and  restricts 
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the  polynomial  coefficients  to  be  non-negative.  L.  J.  Bain 
considers  the  seune  problem  and  hypothesis  testing  in  "Analysis 


for  the  Linear  Failure  Rate  Life  Testing  Distribution" , 

1974  [2U. 

Spline  functions  (segments  of  polynomials)  may  also  be 
used  to  estimate  the  failure  rate  function  or  its  integral. 

The  article  "Spline  Functions  in  Data  Analysis",  1974  [22],  by 
S.  Wold  and  the  program  "FITLOS,  a FORTRAN  Program  for  Fitting 
Low  Order  Polynomial  Splines  by  the  Method  of  Least  Squares” 
by  Patricia  Smith,  1971  I23l,  will  provide  assistance  in  curve 
fitting.  The  use  of  spline  functions  is  appropriate  only 
when  there  is  a large  amount  of  data. 
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3.  APPLICATIONS  TO  ENGINE  MANAGEMENT  PROBLEMS 

The  article  "Classes  of  Distributions  Applicable  in 
Replacement  with  Renewal  Theory  Implications"  by  A.  W,  Marshall 
and  Frank  Proschan,  1972  [8  ] , indicates  the  limitations  of  our 
knowledge  about  optimal  replacement  policies  for  items  with 
life  distributions  in  restricted  classes.  These  limitations 
lead  to  the  need  for  simulation  of  the  engine  management  sys- 
tem to  determine  whether  policies  are  satisfactory.  (The 
simulation  proposed  in  Section  IV  is  one  such  attempt.)  How- 
ever, there  have  been  several  analyses  of  moderately  complex 
machine-repair  and  logistics  systems  that  may  be  of  value  as 
approximations  for  engine  management  planning  and  evaluation. 

The  thesis  "The  Statistical  Analysis  of  Semi-Markov 
Processes  with  Applications  to  Queuing  Problems"  by  T.  L.  Goyal, 
1970  (24  "J,  is  cited  here  because  many  of  the  subsequent  cita- 
tions are  special  cases  of  Semi-Markov  processes.  This  thesis 
(and  similar  articles)  provides  the  statistics  to  estimate 
parameters  of  models  in  the  articles  to  follow. 

The  next  three  citations  are  related  to  replacement  or 
renewal  processes  with  inspection  such  as  the  engine  replace- 
ment process.  J.  M.  Cockerham,  "Maintainability  Analyses  and 
Prediction  Technique",  1967  [25]  , provides  an  analysis  of  an 
"on-off"  process  with  several  modes  of  "off"  depending  on 
whether  an  item  was  removed  from  service  during  operation  or 
at  an  inspection.  The  results  are  limited  to  exponential 


lifetimes.  G.  L.  Lind,  "A  General  Procedure  for  Maintainability 
Prediction",  1968  126],  gives  an  estimation  procedure  for  the 
life  distribution  for  a system  of  components  of  which  some 
fail  in  service  and  some  have  failure  detected  at  the  times 
of  periodic  inspection.  As  in  the  Cockerham  article,  the 
lifetime  distribution  of  components  is  assumed  to  be  exponen- 
tial and  the  system  is  assumed  to  be  parallel-series  in  con- 
struction, i.e.  there  are  spares  for  some  of  the  components. 

The  article,  "Adaptive  Statistical  Procedures  in  Reliability 
and  Maintenance  Problems"  by  J.  L.  Gastwirth  and  J.  H.  Venter, 

1964  [27],  is  representative  of  research  on  inspection  plans. 

The  authors  assume  that  the  lifetime  distribution  is  expontial 
and  derive  an  adaptive  procedure  to  estimate  simultaneously  the 
unknown  parameter  of  the  exponential  distrib«:'cion  and  the 
optimal  inspection  times.  The  adaptive  feature  is  relevant 
for  engine  types  just  going  into  service  for  which  there 
is  little  data.  The  article  is  also  worthwhile  because  it 
indicates  that  inspection  times  as  well  as  maximum  operating 
times  are  subject  to  optimization  as  part  of  maintenance 
policies. 

The  three  articles  cited  in  the  previous  paragraph  have 
application  to  the  engines  themselves  rather  than  the  supply, 
repair,  and  distribution  system.  As  mentioned  repeatedly  in 
Section  1,  any  policies  regarding  engine  management  must  be 
integrated  with  policies  for  operation  of  the  entire  engine 
logistics  system  and  with  Air  Force  flying  hour  predictions.  Articles 
cited  next  present  a broader  view  of  the  engine  management  system. 
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The  article,  "A  Decision  Making  Model  for  Applications 
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of  Queuing  Theory"  by  R.  T.  Nelson,  1970  [28] , is  not  the 
only  article  to  recognize  the  value  of  queuing  analyses  is 
for  decision  making.  However,  it  describes  a queuing  model 
within  a decision  structure  that  could  be  a simplified  repre- 
sentation of  future  engine  management  systems. 

"Generalized  Maintainability  Method  (GMM) " by  R.  M.  Fraser, 
1971  [29],  is  an  extremely  detailed  procedure  that  may  be  applied 
to  the  estimation  of  engine  repair  times . 

Several  branches  of  the  Department  of  Defense  have  sup- 
ported research  on  "multi-echelon"  logistics  systems.  A 
multi-echelon  system  represents  the  engine  repair  system 
which  has  repair  and  inventories  at  bases  and  at  a depot. 

The  final  two  citations  are  the  latest  studies  of  multi- 
echelon systems. 

The  Aerospace  Research  Laboratorj.es  technical  report  by 
H.  P.  Galliher  and  R.  C.  Wilson",  1974  [30],  attempts  to 
justify  the  exponential  assumption  for  the  engine  flying 
hour  distribution.  It  also  includes  simulation  of  several 
variations  on  "METRIC"  inventory  policies.  (METRIC,  Multi- 
Echelon  Technique  for  Recoverable  Inventory  Control,  is  a 
reorder  policy  that  places  an  order  whenever  an  item  is  re- 
moved from  inventory.)  This  report  is  a study  of  the  same 
logistics  system  constituting  part  of  the  engine  management 
system  described  in  Section  I. 
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Last,  but  not  least,  is  the  thesis  "A  Queuing  Theoretic 


Analysis  of  Logistics  Repair  Models  with  Spare  Units"  by 
Francois  bureau,  1974  1 31] . It  represents  the  state  of  analysis 
for  METRIC  type  logistics  systems.  Most  analyses  assume  events 
occur  as  Poisson  processes,  and  this  is  not  necessarily  true. 

It  should  be  indicated  that  even  though  engine  flying  hours 
till  replacement  may  be  approximately  exponentially  distribu- 
ted, this  does  not  mean  that  the  calendar  time  till  replace- 
ment has  an  exponential  distribution.  Remember  that  engine 
life  is  measured  in  flying  hours  and  is  related  to  calendar 
time  only  through  the  flying  hour  program,  whereas  all  other 
activities  in  the  engine  management  system  take  place  in 
calendar  time.  Despite  this  warning,  the  analyses  in  the 
bureau  thesis  provide  approximate  results,  and  the  Poisson 
type  METRIC  logistics  system  analyses  may  be  incorporated 
into  simulation  models  by  the  "control  variates"  technique, 
Fishman  [13 ] , to  sharpen  the  accuracy  of  engine  logistics 
systems  studies. 
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